Inhomogeneous Pairing in Highly Disordered s-wave Superconductors 
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We study a simple model of a two-dimensional s-wave superconductor in the presence of a random 
potential in the regime of large disorder. We first use the Bogoliubov-de Gennes (BdG) approach to 
show that, with increasing disorder the pairing amplitude becomes spatially inhomogeneous, and the 
system cannot be described within conventional approaches for studying disordered superconductors 
which assume a uniform order parameter. In the high disorder regime, we find that the system breaks 
up into superconducting islands (with large pairing amplitude) separated by an insulating sea. We 
show that this inhomogeneity has important implications for the physical properties of this system, 
such as superfluid density and the density of states. We find that a finite spectral gap persists 
in the density of states for all values of disorder and we provide a detailed understanding of this 
remarkable result. We next generalize Anderson's idea of the pairing of exact eigenstates to include 
an inhomogeneous pairing amplitude, and show that it is able to qualitatively capture many of the 
nontrivial features of the full BdG analysis. Finally, we study the transition to a gapped insulating 
state driven by the quantum phase fluctuations about the inhomogeneous superconducting state. 

PACS numbers: 74.20.Mn 74.30. +h 74.20.-z 71.55.Jv 



T3 

S3 
O 

o 



> 
o 

CO 
CN 

O 
O 

03 



i 

S3 
O 

o 



X 



I. INTRODUCTION 

Studies of the interplay between localization and su- 
perconductivity in low dimensions got a boost with the 
possibility of growing amorphous quench condensed films 
of varying thicknesses and degrees of microscopic disor- 
der and the ability to measure their transport properties 
in situ in a controlled manner ^,^]. These experiments 
show a dramatic reduction in T c with increasing disorder 
and eventually a transition to an insulating state above 
a critical disorder strength beyond which resistivity in- 
creases with decreasing T. The data in the vicinity of 
the transition often seems to exhibit scaling behavior, 
suggesting a continuous, disorder-driven superconductor 
(SC) to insulator (I) quantum phase transition at T = 0. 

The physics of these highly disordered films is outside 
the domain of validity of the early theories of dirty su- 
perconductors, due to Anderson || and to Abrikosov and 
Gorkov 0] , which are applicable only in the low disorder 
regime where the mean free path is much longer than 
the inverse Fermi wave- vector. The effect of strong dis- 
order on superconductivity is a challenging theoretical 
problem, as it necessarily involves both interactions and 
disorder §. 

Several different theoretical approaches have been 
taken in the past. First, there are various mean field ap- 
proaches which either extend Anderson's pairing of time- 
reversed exact eigenstates or extend the diagrammatic 
method to high disorder regimes; see, e.g., refs. [p|-pd|. 
In much of the present work, we will also make use of 
mean field theories, which however differ from all previ- 
ous works in a crucial aspect: we will not make any as- 
sumption about the spatial uniformity of the local pairing 
amplitude A. Using the Bogoliubov-de Gennes (BdG) 



approach, as well as a simpler variational treatment us- 
ing exact eigenstates, we will show that outside of the 
weak disorder regime, the spatial inhomogeneity of A 
becomes very important. 

The other point of view, primarily due to Fisher and 
collaborators has been to focus on the universal 

critical properties in the vicinity of the superconductor- 
insulator transition (SIT). These authors have argued 
that fermionic degrees of freedom should be unimportant 
at the transition, which should then be in the same uni- 
versality class as dirty boson problem. As we shall see, 
our results on a simple fermionic model explicitly demon- 
strate how the electrons remain gapped through the tran- 
sition, which is then indeed bosonic in nature. The SC-I 
transition will be shown to be driven by the quantum 
phase fluctuations about the inhomogeneous mean field 
state [p|l. 

We begin by summarizing our main results: 

(1) With increasing disorder, the distribution P(A) of 
the local pairing amplitude A(r) ~ (cj(r)cj(r)) becomes 
very broad, eventually developing considerable weight 
near A s» 0. In contrast, conventional mean-field ap- 
proaches assume a spatially uniform A. 

(2) The spectral gap in the one-particle density of states 
persists even at high disorder in spite of a growing num- 
ber of sites with A(r) w 0. A detailed understanding 
of this surprising effect emerges from a study of the 
spatial variation of A(r) which shows the formation of 
locally superconducting "islands" separated by a non- 
superconducting sea and a very special correlation be- 
tween A(r) and the BdG eigenfunctions. 

(3) We generalize the 'pairing of exact eigenstates' to 
allow for an inhomogeneous pairing amplitude, and find 
that it provides qualitative and quantitative insight into 
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the BdG results. 

(4) There is substantial reduction in the superfluid stiff- 
ness and off-diagonal correlations with increasing disor- 
der. However, the spatial amplitude fluctuations (in re- 
sponse to the random potential) by themselves cannot 
destroy superconductivity. 

(5) We study phase fluctuations about the inhomoge- 
neous SC state described by a quantum XY model whose 
parameters, compressibility and phase stiffness, are ob- 
tained from the BdG mean-field results. A simple anal- 
ysis of this effective model within a self-consistent har- 
monic approximation leads to a transition from the su- 
perconductor to a gapped insulator. 



superfluid stiffness Ds 




FIG. 1. Schematic behavior of superfluid stiffness D s and 
energy gap E gap as a function of temperature T and disorder 
V for the model in Eq. (|). For V = 0, both D s and E gap 
vanish at the critical temperature T c as expected. However, 
for T — the behavior is very unusual with D s vanishing 
at a critical disorder V c but E gap remaining finite and even 
increasing with disorder at large V. 

Schematically our main results on the disorder depen- 
dence of the spectral gap and superfluid stiffness are 
summarized in Fig. |l|. We see that while the super- 
fluid density D s decreases with increasing disorder ul- 
timately vanishing at a critical disorder strength, the en- 
ergy gap always remains finite, and shows an unusual 
non-monotonic behavior: it initially decreases with dis- 
order but remains finite and even increases for large disor- 
der. Note the difference between the finite temperature 
transition in the non-disordered case and the disorder 
driven T = transition. The V = transition at T c is 
driven at weak coupling by the collapse of the gap. In 
contrast the T = transition at V c is driven by a van- 
ishing superfluid stifness even though the gap remains 
finite. 

In ref. 0] we first reported in brief some of these re- 
sults and compared them with earlier Quantum Monte 



Carlo (QMC) studies JlSj of the same model. The present 
paper describes these calculations in detail and extends 
the results to much weaker coupling (the earlier work 
was limited to intermediate coupling with the coherence 
length £ of the order of the interparticle spacing). This 
has been made possible both by technical improvements 
in solving the BdG equations self consistently on larger 
lattices, and by the semi-analytical treatment of the pair- 
ing of exact eigenstates. 

The rest of this paper is organized as follows: In sec- 
tion |l| we describe our model for the disordered SC fol- 
lowed by the inhomogeneous BdG mean field method in 
This section also contains the results of the 
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section 

BdG solution with a detailed discussion of the disorder- 
dependence of various physically interesting quantities, 
such as pairing amplitude, density of states, energy gap, 
order parameter and the superfluid density. In section IV 
we develop the pairing of exact eigenstates theory tak- 
ing into account the inhomogeneity of the pairing ampli- 
tude. Phase fluctuations are discussed in section [v| and 
the phase diagram based on our calculations is described 



in section VI. In section VII we make some predictions 



for STM measurements and discuss some open problems 



and conclude in section VIII 



II. MODEL 

We describe a 2D s-wave SC in the presence of non- 
magnetic impurities by an attractive Hubbard model 
with site-randomness using the Hamiltonian, H = Hq + 
Hint where 



H = -t ^2 ( c L c j<? + h - c -) + y^(Vj - m)"»<7 

<ij>,a i,a 
Hint = -\U\ ^n iT 7i a (1) 



where c ic (ci a ) is the creation (destruction) operator for 
an electron with spin a on a site r$ of a square lattice 
with lattice spacing a — 1, t is the near-neighbor hop- 
ping, \U\ is the pairing interaction, m a — c\ a Ci ai and /j, 
is the chemical potential. The random potential Vi mod- 
els non-magnetic impurities. Vi is an independent ran- 
dom variable at each site r,*, uniformly distributed over 
[—V, V], and V thus controls the strength of the disorder. 

Before proceeding, it may be worthwhile to comment 
on the choice of the Hamiltonian Eq. ([j]). The effects 
of Coulomb repulsion are neglected here in a spirit sim- 
ilar to the Anderson localization problem Jj"6|| . In spite 
of this simplification, Anderson localization has had a 
profound impact on disordered electron systems, and a 
complete understanding of interactions in the presence of 
disorder is still an open problem. Similarly the Hamil- 
tonian we study is a minimal model containing the in- 
terplay of superconductivity and localization. For zero 
disorder V — it describes s-wave superconductivity and 
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for \U\ — it reduces to the (non- interacting) Anderson 
localization problem. We therefore feel that it is very 
important to first understand the physics of this simple 
model before putting in the a dditio nal complication of 
Coulomb effects; (see section VII B for further discus- 
sion). 

We next comment on the choice of parameters. We 
have studied the model ([!]) for a range of parameters: 
0.8 < \U\/t < 8, 0.2 < (n) < 0.875 and a wide range of 
disorder on lattices of sizes up to N = 36 x 36. In our pre- 
vious letter jl4j we reported results mainly for \U\/t = 4 
for two reasons: for computational case, and to compare 
with earlier QMC results. Here we focus on weaker cou- 
pling for which the pure (V — 0) system has a superfluid 
stiffness D s much greater than the spectral gap, the sig- 
nificance of which will become clearer below. Specifically, 
we will show results for \U\/t — 1.5 and (n) = 0.875 on 
systems of typical size 24 x 24. We have taken care to 
work on systems with linear dimension larger than the 
coherence length £ Jl7| . 



III. BOGOLIUBOV DE-GENNES MEAN FIELD 
THEORY 

We use the Bogoliubov de-Gennes (BdG) method Jig} ] 
to analyze the Hamiltonian in Eq. (|lj). This is a mean 
field theory which can treat the spatial variations of the 
pairing amplitude induced by disorder fH^| , which turn 
out to be very important in the large disorder regime 
and lead to many new and unanticipated effects . The 
mean field decomposition of the interaction term gives ex- 
pectation values to the local pairing amplitude and local 
density 

*(r i ) = -\U\(c il c i i), {n la ) = {clc la ) (2) 
and yields an effective quadratic Hamiltonian 

H e g = -t ^2 ( c L c jy + h-c.) + y^(Vi - fii)n ia 

i 

^[A^c^+A*^)^] (3) 



|fn(i"z)| 2 = 1 for each rj. The diagonalization leads to 
the BdG equations 



<IJ>,CT 



where /Jj = n+\U\{rii) /2 incorporates the site-dependent 
Hartree shift in the presence of disorder. Here (rii) — 
J2<j( n i,<?)- We discuss in detail the importance of the 
inhomogeneous Hartree shift in section |v| The effective 
Hamiltonian (Q) can be diagonalized by the canonical 
transformations 



n 



(4) 



Here 7 (7^) is the quasiparticle destruction (creation) 
operator and w„(r;) and u n (rj) satisfy J2n\ u n(^i)\ 2 + 
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A 

-c 



u„(ri) 



(5) 



where the excitation eigenvalues E n > 0. £tt n (rj) = 
-tY.^n^i + 5) + (Vi - fii)u n (ri) where 6 = ±x, ±y, 
and Ku n (ri) = A(rj)u n (rj); and similarly for u n (rj). Ex- 
pressing c i(T ((Her) in terms of quasiparticle operators using 
Eq. (^) the local pairing amplitude and number density 
at T = defined in Eq. (g) are written in terms of the 
eigenvectors of the BdG matrix as 



A(r i ) = \U\J2^n(riK(r i ) 

n 

(n 4 )=2^K(r,)| 2 



(6) 



We now solve the BdG equations, on a finite lattice 
of N sites with periodic boundary conditions, as follows. 
Starting with an initial guess for the local pairing am- 
plitude {A(ri)} and the local chemical potential {/Z^} 
at each site, we numerically determine (using standard 
LAPACK routines) the eigenvalues E n and eigenvectors 
(u n (Yi),v n (Yi)) of the BdG matrix in Eq. (jjj). We then 
compute {A(ri)} and {(Hj)} from Eq. (||). If these val- 
ues differ from the initial choice, the whole process is it- 
erated with a new choice of {A(i"i)} and {(rii)} in the 
BdG matrix until self-consistency is achieved and the 
output is identical to the input for {rii) and A(r^) at 
each site within a certain tolerance. The chemical po- 
tential 11 is determined by l/NJ2i n i = where (n) is 
the given average density of electrons. Note that A(r^) 
and u(rj), can be chosen to be real quantities in the 
absence of a magnetic field. 

The number of iterations necessary to obtain self- 
consistency grows with disorder. We have checked that 
the same solution is obtained for different initial guesses. 
All the results are averaged over 12—15 different real- 
izations of the disorder for a given disorder strength. 



A. Local Pairing Amplitudes and Off Diagonal Long 
Range Order (ODLRO) 

We have compared the ground state energy of the BdG 
solution with that obtained by forcing a uniform pairing 
amplitude, and found that the BdG result was always 
lower. In fact, the difference between the BdG and uni- 
form energies increased with disorder, due to the greater 
variational freedom associated with the inhomogeneous 
BdG solution. 

In Fig. ^ we plot the distribution P(A) of the self con- 
sistent local pairing amplitude A(r^) for several values 
of the disorder V. In the absence of disorder the BDG 
solution has a uniform pairing amplitude Aq ~ 0.153t, 
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identical with the BCS value for V = 0. For low disorder 
V = 0. li, the distribution P(A) has a sharp peak about 
Ao, which justifies the use of a homogeneous mean field 
theory (MFT) (as, e.g., in the usual derivation of An- 
derson's Theorem) for small disorder. With increasing 
disorder V ~ It, the distribution P(A) becomes broad 
and the assumption of a uniform A breaks down. With 
further increase of disorder V ~ 2t, P(A) becomes highly 
skewed with weight building up near A « 0. Qualita- 
tively similar behavior was also found for different choices 
of the attraction U. We found that for the same disorder 
strength V, the fluctuations in A(r.;) are larger for higher 
values of the attraction \U\. 



0.5 



N = 24x24, (n)=0.875, U=-1.5t 




in the non-SC state the off-diagonal correlations decay 
to zero at large distances so Aop = 0. It can be shown 
that Aop — / dAAP(A), i.e., it is the average value of 
the local pairing amplitude. Our calculations show that 
Aop which is identical to Ao in the limit V = 0, is sub- 
stantially reduced by disorder as seen in Fig. ^|. 



N = 24x24, (n)=0.875, U=-1.5t 




0.2 



FIG. 3. The distribution of the number of electrons per 
site n(r) obtained from the BdG solution, for various disorder 
strengths. At low disorder the distribution is sharply peaked 
around the average density (n) — 0.875. P(n) becomes broad 
with increasing V and for large disorder evolves towards a 
bimodal distribution with empty and doubly-occupied sites. 



FIG. 2. Distribution of the local pairing amplitude A(r) 
from BdG calculations for various disorder strength. At 
low disorder the distribution P(A) is sharply peaked around 
Ao ~ 0.15, the pure BCS value for \U\ = 1.5*. P(A) becomes 
broad with increasing V and finally at a very large disorder 
gains significant weight near A ~ 0. 

The distribution of the local pairing amplitude P(A) 
should be contrasted with the distribution of local density 
P(n), which is also inhomogeneous with increasing dis- 
order but very distinct as shown in Fig. ^|. As a function 
of disorder it evolves from being sharply peaked about 
the average (n) at low V towards an almost bimodal dis- 
tribution for large V, with sites being either empty (cor- 
responding to high mountains in the random potential 
topography) or doubly occupied (in the deep valleys of 
the random potential). Later, we will also contrast the 
spatial correlations between the local pairing amplitudes 
and the local densities. 

The off-diagonal long range order (ODLRO) is de- 
fined by the long-distance behavior of the (disorder av- 
eraged) correlation function {c\^c\^CjiCj^) — ► A OP /|[/| 2 
for \Ti — Yj | — * oo. In the SC state Aqp is finite whereas 



B. Single Particle Density of States and Energy Gap 

In Fig. |^ we show the behavior of the single particle 
density of states (DOS) given by 

N M = ZX( r <W«> ~ E ») + + E n)\ (7) 

averaged over disorder. It is seen that with increasing 
disorder the DOS piled up at the gap edge in the pure 
SC is progressively smeared out and states are pushed 
to higher energies. However, the gap in the spectrum 
remains finite. 

The energy gap P gap is obtained directly from the 
BdG calculation as the lowest eigenvalue of the matrix 
in Eq. (JsJ) . From Fig. where we plot the evolution of 
E gap with disorder, we see that not only does the energy 
gap remain finite for all disorder strengths, it even in- 
creases at high disorder! In (much of) the remainder of 
this section and in the following section we will present 
a detailed understanding of this surprising effect. 
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N = 24x24 (n) = 0.875 U = -1.5t 
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FIG. 4. Density of States N(u) for three disorder 
strengths V. With increasing disorder the singular pile-up 
at the gap edge smears out pushing states towards higher en- 
ergies. Surprisingly, the spectral gap remains finite even at 
large V. 



large number of sites with A « at high disorder, one 
might have expected the spectral gap to also collapse. 
However, this expectation is based on an (incorrect) iden- 
tification of the average pairing amplitude, or order pa- 
rameter Aop, with the spectral gap E gap . While the two 
coincide at small disorder strengths, we see from Fig. ^ 
that the two show qualitatively different behavior at high 
disorder. It turns out that important insight into this 
puzzling phenomenon can be obtained by looking at the 
inhomogeneities in A(rj) in real space, as discussed be- 
low. 



C. Formation of Superconducting Islands 

In Fig. |we see the evolution of the spatial distribution 
of the local pairing amplitude for a given realization of 
the random disorder potential with increasing strength of 
the random potential. Though the random potential 
is completely uncorrelated from site to site, the system 
generates, with increasing disorder, spatially correlated 
clusters of sites with large A(r,), or "SC islands", which 
are separated from one another by regions with very small 
A(r,-)Q. The size of the SC islands is the coherence length, 
controlled by the strength of the attraction \U\ and the 
disorder V . 
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FIG. 5. Evolution of -E gap and Aop with disorder. At 
V — both are same as expected; both decrease for small 
disorder, but, the large disorder behavior of them are quite 
different. Also note the non-monotonic behavior of E gap . 

We note that this result is counterintuitive. Given the 
broad distribution of pairing amplitude (Fig. ^|) with a 



FIG. 6. Gray-scale plot for the spatial variation of the 
local pairing amplitude A(r) for a given realization of the 
random disorder potential, which is the same in all the panels 
but with increasing strength of the random potential. Note 
that at large V the system generates "SC islands" (dark re- 
gions) with large pairing amplitude separated by an insulating 
"sea" (white regions) with negligible pairing amplitude. 
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In Fig. we show density n(r^) and A(r^) in a gray 
scale plot for a given realization of the random potential 
at a disorder strength V = 3t. As expected the density 
is large in regions where the random potential is low and 
vice versa on a rather local scale. This is emphasized 
by the density-density correlations being extremely short 
ranged, on the scale of the lattice constant a. The pairing 
amplitude, on the other hand, shows structure i.e the 
formation of SC islands on the scale of the coherence 
length £, which is several lattice spacings. (The coherence 
length of the non-disordered system for this choice of \U\ 
is 6 -10). 




V=3t 



(b) 



5 10 15 30 




A(r) 

FIG. 7. (a) Gray-scale plot of of density n, on the lat- 
tice at strength V = 3i for a given disorder realization, with 
darker regions indicating higher densities, (b) Plot of the 
density-density correlation function for the same disorder re- 
alization n(ri)n(rj) as a function of the distance r = |r< — r 3 - 1 
showing that the correlations decay within a lattice constant, 
(c) Gray-scale plot of of pairing amplitude A(n) on the lattice 
for same V and same realization as in (a), (d) The correlation 
function A(rj)A(r.,) for this disorder realization as a function 
of the distance showing that the correlations persist to dis- 
tances of order several lattice spacings, which is the size of 
the SC islands. 

We next ask: where in space are these "SC islands" 
formed ? The will be very important in our understand- 
ing of the finite energy gap at large disorder. By cor- 
relating the locations of the islands with the underlying 
random potential for many different realizations, one can 
see that large A(r) occurs in regions where \Vi — pn\ is 
small. In the limit of strong disorder, this can be viewed 
as sort of particle-hole mixing in real space. Regions cor- 
responding to deep valleys and to high mountains in the 



potential energy landscape contain fixed number of par- 
ticles per site: two on a valley site or zero on a mountain 
site, and as a result the local pairing amplitude vanishes 
in such regions. 
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FIG. 8. (a) Gray-scale plot of the local pairing ampli- 
tude A(rj) for a given realization of the random potential at 
V = St. (bl-b4) Gray-scale plot of \u n {r)\ 2 + \v n (r)\ 2 for the 
lowest four excitations (n = 1, ... 4); the corresponding eigen- 
values E(l),... -E(4) are also shown. Notice that a particle 
added to (or extracted from) the system has a high proba- 
bility of being found in regions where A(r^) is large. As a 
consequence an energy gap persists even for high disorder. 

To get a better understanding of the finite spectral gap 
-Egap, it is useful to study the eigenfunctions correspond- 
ing to the low energy excitations. In Fig. || we show in 
gray-scale plots the local pairing amplitude A(rj) and 
the |u n (r)| 2 + |w„(r)| 2 for the lowest four excited state 
wave functions, for a given realization of disorder at a 
high value of V — it. We immediately notice the re- 
markable fact (which we have checked for many different 
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realizations) that all the low-lying excitations live on the 
SC islands. Therefore it is no surprise that one ends up 
with a finite pairing gap. 

The obvious next question is: why can one not make 
an even lower energy excitation which lives in the large 
"sea", in between the SC islands, where one would not 
have to pay the cost of the pairing gap? The answer 
is that the random potential makes excitations in the 
"sea" even more energetically unfavorable than those on 
the SC islands. Recall from the preceding subsection 
that the "seas" correspond, roughly speaking, to the high 
mountains and deep valleys in the random potential. It 
is not possible to inject an electron into a deep valley as 
those sites are already doubly occupied, and one needs to 
pay a large (potential) energy cost to extract an electron 
from such sites. Similarly, it is energetically unfavorable 
to create an electron on top of a high mountain in the 
random potential, and it is also not possible to extract 
an electron from such sites as there are none. 

Thus the lowest excitations must correspond to either 
injecting or extracting an electron from regions where 
\Vi — p,\ is small. However, as argued above, these are 
precisely regions with large A(ri), and hence these ex- 
citations have a finite pairing gap. (An understanding 
of the non-monotonic behavior of E gap and its eventual 
increase at large V will come from the analysis of Sec- 
tion 
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An immediate consequence of these ideas is that, while 
the islands may be thought of as locally superconducting, 
the sea separating them can be thought of as "insulating" 
with a large gap determined primarily by the random 
potential. This can be tested by the local density-of- 
states at different regions in the highly disordered regime 
as shown in Fig. [if]. 



E. Superfluid Stiffness 

It is important to understand how disorder, and in 
particular the formation of the inhomogeneous ground 
state, affects the phase rigidity. We calculate the super- 
fluid stiffness D s defined by the induced current response 
to an external magnetic field, given by the usual Kubo 
formula |(]] 

— = (~k x ) - A xx (q x = Q,q y — > 0,iw = 0). (8) 

The first term (— k x ) is the kinetic energy along the x- 
direction and represents the diamagnetic response to an 
external field. (In the continuum the first term would re- 
duce to the total density.) The second term is the param- 
agnetic response given by the (disorder averaged) static 
transverse current-current correlation function: 

^(f x (q,r)f x (- q ,o)) (9) 



1 f 1 ' 1 

A xx (q, iu n ) = — dre 1, 
^ v Jo 



The superfluid stiffness calculated within the BdG 
approximation will be denoted by D®; (the symbol 
D s will be used for the renormalized stiffness to be 
defined later on). Using the BdG transformations 
Eq. (^) the kinetic energy can be written as (— k x ) — 

77 (j2 r n w »i( r ) w »i( r + £)V A straightforward calculation 
of A xx , by expressing the current operators in terms of 
the quasiparticle operators 7 and 7^ , yields at T = the 
result: 



A xx {l 1 2 1 iuo n = Q) = 2t 2 ^ 



1 



(E + E>) 



[v'(2 + x)u(2) + v(2 + x)u'(2)} 
x [u(l + x)v'(l) + v(l)u'(l + x) 
-u(l)v'(l + x) - v(l + x)u'(l)] 
+[u v, v u]. (10) 

Here £ is the unit vector along positive x-direction, 
and we have simplified the notation by using unprimed 
(primed) symbols to denote quantities with subscript nl 
(n2), and n = X,Tj = 2. After disorder averaging, we 
recover translational invariance, so that A xx (r,j, v 3 -,, 0) = 
A xx (Yi — rj-,0) One can then Fourier transform to q to 
obtain A xx (q x — 0,q v ,iu) n = 0), which can be shown to 
go like A+Bqy for small q y . We verify this g y -dependence 
in our numerical results and use it to to take the required 
limit in Eq. (§). 



N = 24x24 (n) = 0.875 U=-1.5t 
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where j£(q) is the paramagnetic current, and w r , 
2-nnT. 



FIG. 9. The superfluid stiffness D° s /-k calculated within 
the BdG theory as a function of disorder. The energy gap is 
also plotted for comparison. Note that, at V = 0, D° 3> E ga , p 
a characteristic of weak-coupling superconductors. However, 
at large disorder one finds D® <C E gs , p suggesting a phase 
fluctuation dominated regime. 
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In Fig. H we show the behavior of the BdG phase stiff- 
ness D®/ir as a function of disorder. The very large re- 
duction of D°, by almost two orders of magnitude, can be 
intuitively understood by the following argument (which 
is also schematically illustrated in Fig. |l0|) . Within mean 
field theory the phases of the order parameter at dif- 
ferent sites are completely aligned in the ground state. 
When an external phase twist 6 is imposed the energy 
of the SC increases leading to a non-zero superfluid stiff- 
ness Dg ~ d 2 E(6)/d8 2 . In a uniform system the exter- 
nal twist is uniformly distributed throughout the system. 
However, in an disordered system where the amplitude 
is highly inhomogeneous the system will distribute the 
phase twists non-uniformly in order to minimize energy 
with most of the twist accommodated in regions where 
the amplitude is small. Thus an inhomogeneous sys- 
tem will be able to greatly reduce its superfluid stiffness 
Ds 0. 



F. Charge Stiffness 

The charge stiffness D° is the strength of S(ui) in the 
optical conductivity <t(uj) and is closely related to D° s . 
It is defined, after analytically continuing A xx to real 
frequency, as pOfl 

D°/tt = (-k x ) - A xx (q = 0,uj^ 0). (11) 

Note the different order of limits compared with the def- 
inition of D s . We found A xx (0;w) ~ co 2 for small oj and 
the resulting value of the charge stiffness is identical with 
the superfluid stiffness: D° = D° for the whole range of 
disorder that we studied. This is to be expected in a 
system with a gap [Q. In fact having established this 
equality, we found it numerically easier to compute D , 
rather than since taking the limit A xx (cu — > 0) is 
simpler than A xx (q y — > 0) on a finite system. 



t t t 



Mean field ground state 




Response to imposed twist 



FIG. 10. (a) Schematic of a disordered SC in which the 
non-uniform amplitude results in the formation of SC islands. 
Within mean field theory the phase in the ground state is 
spatially uniform even though the amplitude is not (larger 
amplitudes is represented by larger arrows), (b) Schematic 
illustration of the response to an externally applied phase 
twist (of 7r/2) indicated by the fat arrows. The system has 
a non-uniform response, with larger phase twists in regions 
where the amplitude is small. This results in a smaller stiff- 
ness D° ~ d 2 E(6)/dd 2 compared to the case of a uniformly 
distributed phase twist. 

We emphasize that despite this dramatic reduction in 
Dg, the superfluid stiffness continues to remain non-zero 
within the BdG approximation. In other words, the fluc- 
tuations in the pairing amplitude alone are unable to 
drive the system into an insulator. In order to describe 
the SC-Insulator transition it is therefore essential to take 
into account phase fluctuations as discussed in Section |V| 
below. 



IV. PAIRING OF EXACT EIGENSTATES 

Although the BdG analysis described in the previous 
section led to various striking results, and considerable 
physical insight, there were a few points which could not 
be adequately addressed: (1) We could not study the 
weak-coupling limit \U\/t -C 1, since the exponentially 
large coherence length £ leads to severe finite size effects 
in the numerical calculations. (2) Although the existence 
of the gap at large disorder could be understood, we did 
not get any insight into its non-monotonic dependence 
on disorder. 

In order to address these issues, and to gain a deeper 
understanding of the BdG results, we now generalize An- 
derson's original idea of pairing the time-reversed exact 
eigenstates of the disordered, noranteracting system [||, 
in a manner that that allows the local pairing amplitude 
to become spatially inhomogeneous. We will show that 
this generalization permits us to recover most, but not 
all, of the qualitative features of the BdG results. This 
analysis also has the virtue of leading to much simpler 
equations from which one can gain qualitative insights in 
the weak and strong disorder limit. 

We begin with the noninteracting disordered problem 
defined by the quadratic Hamiltonian Ho of eq. (|l|) . The 
corresponding eigenvalue problem is, in principle, solu- 
ble: 7io\4> a ) = £a\<fia), where a labels the exact eigen- 
states of H.Q. Next, following Anderson let us imagine 
pairing electrons in time reversed eigenstates a, f and 
a, J.. The analog of the "reduced BCS" Hamiltonian in 
this basis is then given by 

W = ^t, a J aa c a(T - lUl^M^pcl^c^cp! (12) 

a, a a,0 

where the matrix M Q .^ is defined by 

M a ^ = ^|^ a (r i )| 2 |^(r i )| 2 - (13) 
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Here £ Q = (e a — jx) is measured relative to the aver- 
age Hartree shifted jj, which fixes the electronic density. 
(We will return to the question of average versus site- 
dependent Hartree shifts later in this section.) A BCS- 
like mean field decomposition of the Hamiltonian jl^ ) 
leads to the T = gap equation 



A a = \U\J2 M *, 



A/1 
'2Ep 



(14) 



where E a — v/£l + A^ . Further, the chemical potential 
jl is determined by the number equation 



En 



(15) 



Our formulation generalizes Anderson's original anal- 
ysis by retaining the full M a ^, and it is the structure 
of this matrix which will permit us to access the large 
disorder regime with highly inhomogeneous pairing. 



A. Non-monotonic behavior of energy gap as a 
function of disorder 



nearby in space are far in energy and vice-versa. Next 
we identify ^ r |0 Q (ri)| 4 as the participation ratio for 
the (normalized) state <fi a (ri), which in turn is given by 
Cioc (oQ where Cioc(a) is the localization length for that 
state 0. 

Thus for large disorder we solve the gap equation (|l4| ) 
with the kernel M a ss ^Q./3Cioc ( a )- We nn( i that f° r 
states a with energies far from the chemical potential, 
the solution is A Q = 0, i.e., these states are unaffected 
by pairing. On the other hand, for states with small £ Q 
we find E a ~ |t/|/[2(^ c (a)]. One thus obtains a gap: 



E — 



Sloe 



(16) 



in the high disorder limit, where (io C is the localization 
length at the chemical potential. 

The diagonal approximation becomes exact in the ex- 
treme site localized limit (V — > oo). In this case, the role 
of the exact eigenstate label a is played by the site i\ 
at which the state is localized. It is easy to show that 
all states for which £ ri < \U\/2 have finite pairing ampli- 
tude and a spectral gap of |U|/2, which is a well known 
result 



We now solve the equations obtained above from the 
pairing of exact eigenstates. A qualitative analysis of the 
large and small disorder limits will be given and com- 
pared with a full numerical solution, and all of these re- 
sults will be compared with the BdG results of the pre- 
vious Section. 

Let us begin with the low disorder regime. For a finite 
system in 2D, or an infinite system in 3D, the eigenstates 
4>a(^i)'s are extended on the scale of system. We thus 
find M Qj/ 3 w 1/N, independent of a and /?, which we call 
the "uniform approximation" for M. This is this limit 
in which Anderson's theorem applies, the gap equation 
takes the simple BCS form, and one obtains a (spatially) 
uniform A. 

The behavior of E gap with this simplified picture is cal- 
culated and shown as "uniform approximation" in Fig. [l^ 
for low V. The decrease of -Egap with increasing V in 
this regime can be traced primarily to a simple density- 
of-states effect in the BCS result for the gap. For the 
nearest-neighbor dispersion in 2D and the filling chosen, 
one can easily show that the average DOS at the chem- 
ical potential, N(£ = 0), decreases with increasing V in 
the weak disorder limit. 

In the high disorder regime, on the other hand, the 
eigenstates of the non-interacting problem are strongly 
localized and different states have a very small spatial 
overlap. We can therefore make a "diagonal approxima- 
tion" for for the M-matrix: M a _p fa S a _/3j2 r - I0a( r i)| 4 - 
We have numerically checked that the diagonal ele- 
ments of M are indeed the largest elements as shown 
in Fig. (p~T]) . Moreover, the off-diagonal elements are not 
important in the gap equation (J14|) , as states which are 



FIG. 11. Grey-scale plot of the matrix elements of M at p 
at large disorder V = 6t for a 30 x 30 non-interacting system. 
The x and y axis are the a and j3 indices respectively. Note 
that diagonal matrix elements are the largest. 

In Fig. [l^ we compare the large disorder asymptotic re- 
sult from Eq. ([l6]) , the "diagonal approximation" , with a 
full numerical solution of equations (^) and ([if]) of the 
method of pairing of exact eigenstates (where we self- 
consistently determined A a 's for all a's and /}). Finally 
we also show in Fig. [l^ the BdG solution, with a uni- 
form Hartree shift, which is in excellent agreement with 







the results from pairing of exact eigenstates. (Similar 
agreement is also found for all the other quantities like 
the P(A), A OP , N(u) and D° s /n as a function of V) 

To summarize: we now have a complete understanding 
of the non-monotonic dependence of the spectral gap on 
disorder. The weak disorder asymptotic shows that the 
initial drop is a simple density-of-states effect. On the 
other hand the increase of the gap in the strong disorder 
limit comes from the decrease in the localization length 
Cioc as seen from Eq. ([To]). 

It is important to emphasize that while the numeri- 
cal comparisons in Fig. [l2| are for a moderate value of 
\U\ = 1.5t, the method of pairing of exact eigenstates 
should work best in the weak-coupling limit, where \U\ is 
the smallest energy scale in the problem, and hence the 
noninteracting problem is diagonalized first. The analyt- 
ical approximations in the small and large disorder limits 
given above are thus valid even for \U\/t <C 1 where we 
cannot do reliable numerical calculations. 
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FIG. 12. Upper Panel : Comparison of the energy gap 
-Egap as a function of disorder obtained by the generalized 
pairing of exact eigenstates method and the BdG approach. 
Both methods are implemented with only an average Hartree 
shift. Also shown are two asymptotic regimes for the gap. 
The decrease of E gap at low V because of the decrease of the 
DOS with increasing V is described within a "uniform approx- 
imation" (see text) and the increase of -E gap at large disorder 
due to strong localization effects on the single particle wave 
functions is described within the "diagonal approximation" 
(see text) . Lower Panel : Comparison of the energy gap -E gap 
as a function of disorder calculated within the BdG approach 
with an average and a site-dependent Hartree shift. While 
the two results are qualitatively similar, there are quantita- 
tive differences. 



B. Correlation between SC islands and excitations 

The pairing of exact eigenstates formulation also gives 
analytical insight into the large spatial overlap of the 
low lying excited states with the SC islands in the large 
disorder regime, which was observed and discussed at 
length in the previous section. 

One can show quite generally that the pairing am- 
plitude in real space A(r) is related to A(a) through 
A Q = J2 Ti A(r i )|0 Q (r l )| 2 . The gap equation Q can 
then be rewritten as 



A(rO 



\U\ 
2 



E 



= l0a(rO| s 



(17) 



We now specialize to the large disorder regime and use 
the solution of the gap equation within the "diagonal ap- 
proximation" in the preceding subsection to note that 
the only a'a which contribute to the sum are those with 
£o ~ 0, since otherwise A Q = 0. The above equation then 
simplifies to A(r^) w \U\ J2' a |^a(i"i)| 2 /2 with the sum re- 
stricted to states near the chemical potential. This im- 
mediately shows the strong correlation between the spa- 
tial structures of the regions of A(r^), the SC islands, and 
that of the eigenstates <5f> a ( r i) of the non- interacting prob- 
lem, which are are precisely the low-lying excitations. 



C. Importance of site-dependent Hartree shift 

Having seen the great success of the pairing of exact 
eigenstates in reproducing the results of the BdG anal- 
ysis, we finally turn to the one important feature of the 
BdG analysis that is not captured in the present frame- 
work. We saw in Eq. (||) that the BdG equations incor- 
porate site-dependent Hartree shifts, while the method 
of exact eigenstates did not. We now discuss what the 
effects of inhomogeneous Hartree terms are and why such 
terms are not easy to deal with in the exact eigenstates 
formalism. We note that we are not aware of any previous 
work that has looked at the effects of such inhomogeneous 
Hartree shifts. 

First, inclusion of the site-dependent Hartree terms 
lead to a quantitative change in the behavior of £gap as 
a function of V as seen from the lower panel of Fig. |lj. 
The calculation with the uniform (average) Hartree shift 
shows qualitatively similar non-monotonic behavior with 
a minimum in the gap at a slightly smaller value of V. 

A much more dramatic effect can be seen in the den- 
sity of states plotted in Fig. (|l3|). The calculation with an 
average Hartree shift has a BCS-like pile-up in the DOS 
at the gap edge, while the result with the site-dependent 
shifts shows that this pile-up is completely smeared out 
and states are pushed out to the band tails. The occur- 
rence of the DOS peak within the theory of pairing of 
exact eigenstates (with homogeneous Hartree shift) has 
the same origin as in BCS theory. However, all ener- 
gies are measured with respect to chemical potential, and 
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when the inhomogeneity in the Hartree shift of chemical 
potential is taken into account, it acts like a random per- 
turbation which breaks the degeneracy of states near gap 
edge. 

It would have been nice to incorporate a site-dependent 
Hartree shift in the theory of pairing of exact eigen- 
states. However, in this case the "normal state" Hamil- 
tonian whose exact eigenstates one would have to solve 

for would be: ^Normal = -tJ2<ij>,a( c ia c 3°- + ^- c -) + 

J2i — L 1 — \U\(rii)/2) n ia . However, one then loses 
much of the simplicity of the exact eigenstates formal- 
ism since ^Normal is itself an interacting problem, which 
needs to be solved self-consistently. Further, there are 
some problems (which we will not discuss here) associ- 
ated with treating U at the Hartree level alone, before 
incorporating the pairing physics, in the large disorder 
regime p2| . 

In conclusion, while the generalized pairing of exact 
eigenstates is able to give much insight into the behavior 
of the spectral gap and pairing amplitudes, and gives 
qualitative information about the weak coupling limit, 
the BdG method with site-dependent Hartree shifts is 
the best scheme for quantitative results. 
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FIG. 13. The density of states obtained by using the pair- 
ing of exact eigenstates method using only the average Hartree 
shift at a disorder strength V — 2t in the insulator. Note the 
unphysical pile up of states just beyond the gap. This artifact 
is not present in the BdG density of states shown in Fig. (^) 
which includes the local Hartree shift also self consistently. 



V. QUANTUM PHASE FLUCTUATIONS 

As described at the end of Section III, the BdG analysis 
leads to a large suppression of the superfluid stiffness, 



but the disorder induced amplitude inhomogeneity is not 
sufficient to drive D s to zero. In order to understand the 
transition to an insulating state, we must focus on the 
phase degrees of freedom which are ignored (or frozen) 
in the mean field description used thus far. We thus use 
the 2D quantum XY action in imaginary time: 



Sg = ^l dr^ 



2 r 



V dr 



4 



/ dTV(l-co S [0(r,r)-0(r + J,T)]) (18) 

JO ,. A 



to describe the dynamics of the phase variables 9(r, r) de- 
fined on a coarse-grained square lattice of lattice spacing 

We can motivate the use of such a model in both the 
weak disorder limits and therefore use it for all disorder 
strengths. At weak disorder, one can follow the deriva- 
tion of ref. H to derive an effective action for the phase 
variables in a disordered systems, and then coarse grain 
to the scale of £ using the method of ref. p3] to obtain 
the above action. This coarse graining shows that the 
coefficient of the time derivative term is £ 2 k in 2D where 
K = dn/dfi is the static, long- wavelength compressibility 
calculated at the mean-field level, and the coefficient of 
the cosine term is the mean-field phase stiffness D®. 

In the opposite high disorder limit one can view Jig) as 
simply describing a Josephson junction array of the SC 
islands embedded in an insulating sea as seen in Fig. |(]. 
In this case, the first term represents the charging energy 
of the islands and the second term the Josephson cou- 
pling between islands. Further we are making the crude 
approximation of ignoring the random variations of the 
charging and coupling energies in this random system, 
and simply using the mean field values obtained from 
the BdG analysis. We also ignore the disorder depen- 
dence of the coherence length £, and for simplicity use 
its V = value £ - 

The nonlinearities in the cosine term lead to a renor- 
malization of the superfluid stiffness which we calcu- 
late within the self-consistent harmonic approximation 
(SCHA) |Q. This is done by choosing a trial Gaussian 
action 



So 



fd9{r, t) 



V dr 

+ 1T I" drJ2(e^r)-e(r + S,r)r (19) 

8 J ° r,6 

where the renormalized stiffness D s is determined by us- 
ing the following variational principle 



Fg < F + (Sg - So) . 



(20) 



Here Fg and F are the free energies corresponding to 
the actions Sg and So respectively, and the expectation 
value (. . .)o is determined using the Gaussian Sq. 
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Minimizing the right hand side of eq. (|20|), we obtain 
the result 

D s =£> s °exp(-(C.) /2). (21) 



Here (0f,-)o is the mean square fluctuation of the near- 
neighbor phase difference, and is given by 

NX 



E 

Q 



(22) 



Here sq = 2 [2 — cos((2:r) — cos(Q 1/ )], and the momentum 
sum is restricted to Qi < it. Defining the renormalization 
factor X = D s /D°, and 



D o <2 I jvZ, Q J 



(23) 



which depends upon all of the bare parameters of the 
theory, we find from eqs. (||) and d|) that X is given 
by the solution of 



X = cxp I — y/a/X 



(24) 



We solve this transcendental equation to determine the 
renormalized D s a function of V, using as input for a the 
BdG results for the bare stiffness D° and compressibil- 
ity k for each value of V. The BdG compressibility is 
plotted in Fig. |l4|(a). We do not have a simple physical 
for the small maximum in k at low disorder, which is a 
parameter dependent feature absent for larger values of 
\U\. However, our results for the renormalization of D s 
are completely insensitive to the presence or absence of 
this non-monotonicity. 

The renormalized D s obtained from the SCHA is plot- 
ted in Fig. |lj(b) as the full line. Quantum phase fluc- 
tuations are found to lower the stiffness and beyond a 
certain critical disorder drive it to zero, even though 
the bare (BdG) stiffness was always nonvanishing. Thus 
the SCHA gives a transition to a non-superconducting 
state, even though it is unreliable in vicinity of the tran- 
sition. In particular, it predicts a first order transition at 
ctcrit = 4exp(— 2) with a jump discontinuity of exp(— 2) in 
the value of X. We believe that the order of the transition 
is an artifact of the approximation, although the critical 
disorder obtained from such a calculation is in reason- 
able agreement with quantum Monte Carlo results jig] 
for parameter values (\U\/t = 4) for which a comparison 
can be made JTJ] . 

We next argue that quantum phase fluctuations do not 
have a significant effect on the electronic excitation spec- 
trum. This is because the spectral gap at large disorder 
arises from low energy excitations that live on a SC is- 
land, which is relatively unaffected by phase fluctuations. 
On the other hand, as we have seen above, these fluctua- 
tions have a profound effect on suppressing the coherence 
between SC islands. Thus the non-superconductor in this 



model continues to have a finite spectral gap for one- 
electron excitations even after the effects of phase fluc- 
tuations arc included. Thus it is an insulating state. Fi- 
nally the absence of low lying electronic excitations near 
the transition implies that the quantum phase transition 
in this electronic model is in the superfluid-Bose insulator 
universality class fl2|. 
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FIG. 14. (a) Left panel: The compressibility k — dn/dfi 
as a function of disorder V. (b) Right Panel: Evolution of 
superfluid stiffness D 3 /n upon including the quantum phase 
fluctuations. The bare BdG stiffness D° is plotted as symbols 
with a dashed line through them, while the renormalized stiff- 
ness D s /n is shown by the full line. D s vanishes at V c = 1.75t 
beyond which the system is insulating. 



VI. PHASE DIAGRAM 

In this section we discuss the T — phase diagram 
for the disordered, attractive Hubbard model in the 
(|E7|/i,V/i)-plane. It is known @ that, on the \U\ = 
axis, for all values of disorder F / 0, one has an An- 
derson insulator with gapless excitations in 2D. On the 
V = axis one simply has a crossover as a function of 
\U\/t from a BCS superconductor to a condensate of in- 
teracting (hard core) bosons |25| ]. 

The four symbols marked in Fig. |l5| are the result of 
BdG analysis supplemented by the simple phase fluctua- 
tion analysis described above. Despite the simplifying ap- 
proximations involved, and the lack of a detailed study of 
finite size effects, we nevertheless believe that our results 
do give a reasonable qualitative idea about the critical 
disorder V C (U) separating the SC phase from an insula- 
tor with a gap in its single-particle excitation spectrum. 
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Further our estimated V c at \U\Jt = 4 is in reasonable 
agreement Jl4| with quantum Monte Carlo results (ljj . 
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FIG. 15. Schematic phase diagram at T = of the disor- 
dered, attractive Hubbard model in the disorder V— attraction 
| U plane. The entire y-axis (\U\/t — 0) corresponds to an An- 
derson insulator with gapless excitations. At finite \U\/t there 
are two phases- a superconducting phase at low disorder and 
a gapped insulating phase at high disorder. Thus U is a sin- 
gular perturbation in that the smallest \U\ induces a gap. The 
symbols denote the critical disorder V C (U), separating the su- 
perconducting and the gapped insulating phases, estimated 
from the calculations described in the text. We argue against 
possibilities (b) and (c) for the form of the phase boundary 
in the \U\ — > limit, and suggest that V C (U —* 0) approaches 
a finite value of order unity, as shown schematically by curve 
(a). We find no evidence for a gapless Fermi insulator phase. 



In principle, there are three possibilities for the con- 
tinuation of the V C (U) phase boundary as \U\/t — > 0, a 
limit which we cannot address numerically. As shown in 
Fig. [l5| these are: (a) V C (U — ► 0) is a finite number of 
order unity; or (b) V C (U — ► 0) diverges to infinity; or (c) 
V C (U — ► 0) vanishes. We will now argue against (b) and 
(c), suggesting that (a) is in fact the correct result. 

First we examine possibility (b) by looking at the case 
of a fixed small \U\/t with V — > oo. ^From the large 
disorder asymptotics of the preceding Section (within the 
"diagonal approximation" for the matrix M) we found 
that one obtains SC islands whose size is the localization 
length. Thus the effective coherence length is determined 
by Cioc, he., the disorder and not by the weak coupling. 
Since this length scale becomes very small for large V, we 
expect phase fluctuations to destroy the long range phase 
coherence between the small SC islands. Thus we find it 
very hard for SC to persist out to very large disorder as 



required by the possibility (b). 

Next consider possibility (c) by studying the case of a 
fixed, small V taking the limit \U\/t <C 1. Here one can 
just use the standard theory of dirty superconductors. 
The pure (V = 0) coherence length £o is exponentially 
large in \U\/t, and even if the coherence length in the 
disordered problem is given by £ ~ \/£o^ £ neverthe- 
less grows as \U\/t is reduced. With a growing coherence 
length, both amplitude and phase fluctuations are sup- 
pressed, and we cannot see how SC can be destroyed as 
required by the possibility (c). 

There have been suggestions from QMC studies 
of two insulating phases: a gapless "Fermi" insulator at 
small \U\ and a gapped "Bose" insulator at large \U\ for 
the model in Eq. (|l|). It is possible that a vanishing gap 
may have been observed because of the finite temperature 
in the simulations. We see absolutely no evidence for a 
"Fermi" insulator in this model, away from the \U\ = 
line, and we have presented strong numerical evidence 
and arguments for a finite gap in the non-SC state for 
any \U\ > 0. 

In the \U\/t 1 our Hamiltonian maps on to the 
problem of hard core interacting bosons, with an effec- 
tive hopping ti, ose ~ t 2 /\U\, in a random potential. For 
this problem one expects 14(|?7| — > oo) ~ t 2 /\U\, which 
gives us an understanding of the decrease in V c with \U\. 
Further, in this limit the insulating phase is precisely the 
Bose glass phase |12] . 



VII. FUTURE DIRECTIONS 



A. Prediction for STM measurement 



As discussed in Sec. ( III C), we find that the disordered 
SC consists of "SC islands" with significant pairing am- 
plitude which are separated from the insulating sea with 
nearly zero pairing amplitude. Further we predict that 
the local energy gap in the SC islands (upper panel) 
should be smaller than the local gap in the insulating 
sea (lower panel) as seen by the nature of the local den- 
sity of states (LDOS) in Fig. (jig). The insulating nature 
of the sea is brought out by the absence of any build up 
of weight in the LDOS near the gap edge. It should be 
possible to measure the LDOS using an STM probe as 
has already been demonstrated in the context of mag- 
netic impurities in s-wave superconductors [ ^7f and with 
non-magnetic impurities in the high T c d-wave supercon- 
ductors EiS. 



B. Coulomb Effects 

Coulomb interactions have been typically included as 
an effective //* term in the coupling p). However, we ex- 
pect that Coulomb interactions coupled with an inhomo- 
geneous pairing amplitude should produce qualitatively 
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novel effects. The effective attraction between electrons 
could become inhomogeneous producing regions where 
the pairing amplitude is suppressed with locally gapless 
excitations. Another possibility is that in the presence 
of repulsion between electrons, local moments could be 
formed |30| in the disordered SC which could act like 
magnetic impurities and be pair-breaking. This is an im- 
portant direction for further study. 



N = 30x30, (n)=0.B75, U=-2t, V=2.5t 




FIG. 16. (a) Upper panel: The local density of states 
(LDOS) at sites where the pairing amplitude A is large. These 
regions correspond to the "SC-islands" which have a small lo- 
cal superconducting gap and a coherence peak at the gap 
edge, (b) Lower Panel: LDOS at sites with A ~ 0. These 
regions correspond to the insulating sea showing a larger spec- 
tral gap. Also the coherence peak features at the gap edge 
are washed out as a result these regions show a pseudogap 
behavior. 



whereas the transition in the granular films of category 
(b) is driven by the loss of phase coherence. In granular 
films superconductivity sets in in two steps: with decreas- 
ing temperature first the individual grains become super- 
conducting but the Josephson coupling between grains is 
still weak so the film is still resistive. As the temper- 
ature is lowered the Josephson coupling is able to lock 
the phases of the individual grains into a globally phase 
coherent state at the transition temperature T c below 
which the film resistance is zero. 

We have shown from our calculations that inhomoge- 
neous ramified structures are generated even in models 
which are "homogeneously" disordered. 

We find that a 2D homogeneously disordered SC shows 
a highly inhomogeneous pairing amplitude. At higher 
disorder, we obtain superconducting "islands" or grains 
which are coupled by Josephson coupling. Thus our the- 
oretical work, as described in this paper, challenges the 
picture of two different paradigms to describe the SC- 
insulator transition - an amplitude-driven mechanism for 
the nominally 'homogeneous films' and a phase-driven 
mechanism for the granular films. We find instead that 
both amplitude and phase inhomogeneities play an im- 
portant role. With increasing disorder first the amplitude 
inhomogeneity is important in generating weak links cor- 
responding to small pairing amplitude. These weak links 
then become the active sites where Josephson coupling 
between the SC islands is weakened and leads to a de- 
coupling of the phases between islands and hence a SC-I 
quantum phase transition. 

Our calculations support the conjecture of Fisher et. 
al. (REFERNCE) that the SIT is in the universality 
class of disordered bosons, since we find that once the 
pairing amplitude is allowed to become inhomogeneous, 
there are no low lying fermionic excitations in the disor- 
dered s-wave superconductor. The transition to a gapped 
insulator is driven by the vanishing of the superfluid stiff- 
ness due to quantum phase fluctuations about the inho- 
mogeneous paired state. 
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VIII. CONCLUSIONS 

Depending on the material of the film and the type 
of substrate it is experimentally possible to grow two 
types of films- (a) nominally homogeneously disordered 
films pi)] , and (b) granular films [|32|]33| 1. It is believed 
that the nature of the SC-insulator transition in these 
two types of films (either as a function of temperature 
for a given film, or as a function of film thickness at 
the extrapolated T = transition for a sequence of film 
thicknesses), is quite distinct. The transition in the so- 
called 'homogeneous' films of category (a) is believed to 
be driven by the collapse of the superconducting ampli- 
tude as a function of increasing temperature or disorder, 
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